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Abstract 

In this paper, we study the efficient numerical integration of functions with sharp 
gradients and cusps. An adaptive integration algorithm is presented that systematically 
improves the accuracy of the integration of a set of functions. The algorithm is based 
on a divide and conquer strategy and is independent of the location of the sharp gradi- 
ent or cusp. The error analysis reveals that for a C° function (derivative-discontinuity 
at a point), a rate of convergence of n + 1 is obtained in W 1 . Two applications of the 
adaptive integration scheme are studied. First, we use the adaptive quadratures for the 
integration of the regularized Heaviside function — a strongly localized function that is 
used for modeling sharp gradients. Then, the adaptive quadratures are employed in 
the enriched finite element solution of the all-electron Coulomb problem in crystalline 
diamond. The source term and enrichment functions of this problem have sharp gradi- 
ents and cusps at the nuclei. We show that the optimal rate of convergence is obtained 
with only a marginal increase in the number of integration points with respect to the 
pure finite element solution with the same number of elements. The adaptive integra- 
tion scheme is simple, robust, and directly applicable to any generalized finite element 
method employing enrichments with sharp local variations or cusps in n-dimensional 
parallelepiped elements. 



1 Introduction 

Functions with sharp gradients appear in the solution of problems with localization and 
cohesive process zones [1-3], shear bands [4,5], thermal gradients [6-8], convection- and 
advection-diffusion problems [9-11], and in electronic structure calculations and ab initio 
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materials modeling [12-14]. In the partition-of-unity [15] solution of these problems, such 
sharp gradients and cusps are efficiently resolved by incorporating enrichment functions that 
resemble the solution locally. As a result, efficient numerical integration of the basis functions 
and their gradients to form the system matrices becomes computationally demanding since 
one has to deal with strongly localized functions, instead of polynomial integrands. In 
many applications, the enrichment functions are the solution of local problems and known 
only numerically. Evaluation of such integrands can be extremely time-consuming, which 
points to the need for an efficient integration scheme. In this paper, we present an adaptive 
scheme for the integration of functions with sharp gradients and cusps. Adaptive integration 
algorithms accumulate integration points in regions with higher errors, and use fewer points 
where the integrand is smooth. Also, we are interested in domains that can be prescribed as 
a collection of hyperparallelepipeds, for example, parallelograms and parallelepipeds in two 
and three dimensions, respectively. 

Adaptive integration schemes are normally recursive in nature and have a few common 
ingredients [16]: a quadrature rule that can be applied to the integration domain to provide 
a local estimate of the integration; a procedure to estimate the local integration error; a 
strategy to partition the integration domain into smaller divisions of the same shape; and 
a stopping criterion. The algorithm of Gander and Gautschi [17] over the interval and that 
of Berntsen et al. [16, 18] for a collection of triangles and tetrahedra are examples. Genz 
and Cools [19] proposed an algorithm for a vector-valued function over a combination of 
n-dimensional simplices. At each step, a subset of the simplices with the highest errors 
are selected and subdivided, and quadratures over the subdivisions are used to update the 
integral. The local error in Reference [19] is estimated by applying null quadrature rules — 
quadratures that (incorrectly) integrate to zero all polynomials up to degree rf, and fail to do 
so for at least one polynomial of a higher degree d+ 1 [20, 21]. Herein, we use the difference 
in the integration of a function with two different tensor-product quadrature rules as the 
local error estimate, and proceed to subdivide a cell until the absolute error of integration 
of all the functions over each individual cell falls below a prescribed tolerance. 

In this paper, we consider two types of local features. First, we focus on the regularized 
Heaviside function: this function has a sharp gradient, and by shrinking the size of the 
process zone, becomes strongly localized. Tornberg [22] and Patzak and Jirasek [1] proposed 
regularized forms of the Heaviside function that smear the strong discontinuity over a short 
distance on which the discontinuous function is approximated by a localized function. Oh et 
al. [23] introduced several smooth-piecewise-polynomial regularized discontinuous functions 
that can be used in the extended finite element method. Benvenuti et al. [2] presented a 
method to integrate the regularized Heaviside function by the integration of the equivalent 
smooth functions. This is an extension of the method of Ventura [24] for the integration 
of discontinuous functions, which is however restricted to elements with constant Jacobian 
of the transformation. A more general treatment for the integration of arbitrary classes of 
functions, which is also adopted in this work, is adaptive integration using a posteriori error 
estimates: the integration points are concentrated close to the region where sharp gradients 
appear, and fewer points are used elsewhere [2]. 

Next, we consider functions with cusps. We pick the Coulomb problem in crystalline 
diamond and apply an enriched finite element (EFE) approach: Poisson's equation is solved 
with the electronic charge density as the source term, and the isolated atom solutions as the 
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enrichment functions [14]. These functions have sharp gradients in the region close to the 
atomic sites, and have cusps at the nuclei. In one dimension, Mobius transformation has 
been used for the integration of functions with a peak at or near a boundary [25-27]. When 
Mobius transformation is applied to a standard quadrature over the interval, integration 
points are attracted toward the sharp gradient and more accurate results are obtained. 
However, this technique cannot be applied in higher dimensions. A well-known remedy in 
such cases is adaptive integration. Our adaptive integration scheme recursively performs a 
uniform refinement of the parent cell until the prescribed error tolerance over each subcell is 
met. The number of integration points in the adaptive quadratures is only marginally more 
than that obtained for the finite element solution with the same number of elements (and 
much lower accuracy) , which points to the fact that numerical integration is not a bottleneck 
for the EFE solution of the problem. 

The structure of this paper follows. In Section 2, the adaptive integration scheme is 
introduced, followed by an error analysis in Section 3. It is shown that the algorithm is 
efficient for integrands with a cusp, such as C° functions of the form /(x) = 1 — r or 
/(x) = exp(— ar). The adaptive integration algorithm is used for the integration of the 
regularized Heaviside function in Section 4.1. In Section 4.2, an enriched finite element 
approach is applied to solve the all-electron Coulomb problem in crystalline diamond (charge 
density has a cusp at the nuclei). We close with a few concluding remarks in Section 5. 

2 Adaptive Integration Scheme 

Our quadrature construction algorithm is customized to meet the integration demands of 
high accuracy over parallelepipeds. A tensor-product quadrature based on a one-dimensional 
Gauss rule with five points in each direction is used to evaluate the local integrals. A tensor- 
product quadrature with eight points in each direction is used to compute the reference 
integral, and provides an estimate of the local integration error. The reason for choosing 
eight points in each direction is that a quadrature with six points might not be accurate 
enough to detect the integration error (due to the sharp gradients of the integrands); and a 
seven-point quadrature may suffer from the same odd-even defect of the quadratures as does 
a five-point quadrature. If the absolute error of integration is greater than the prescribed 
tolerance, the integration domain is uniformly divided into eight cells (in three dimensions) 
and the adaptive integration is performed over each cell recursively. This process is started 
with all the functions that need to be integrated, and at each step only those functions whose 
integration error is greater than the tolerance are passed to the next level, until all functions 
are integrated to the specified tolerance. Using the relative error as the stopping criterion 
may cause numerical difficulties in parts of the domain where the integral of at least one of the 
functions is close to zero. This problem was also reported in Reference [17], and was resolved 
by additionally checking the local estimate against the global estimate of the integration: if 
the contribution of the partition was small enough, further subdivision was circumvented. 
Such a problem is avoided by using the absolute error as the measure of accuracy. Note 
that the error tolerance in our quadrature construction scheme is only a stopping criterion, 
and may not represent the exact integration error. However, it provides a systematic means 
to improve the accuracy of the numerical integration until stable solutions are obtained. In 
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a specific application, the appropriate tolerance depends on the overall accuracy required 
and the mesh resolution: as the mesh is refined, the tolerance (error per element) should 
be decreased proportionally to attain the same overall accuracy. Furthermore, if material 
parameters that appear in the weak form integrals have strong spatial variations, they can 
be included in the integrands fa (see Algorithm below) that are evaluated to construct the 
quadrature rule within the element. The following pseudo-code explains our quadrature 
construction scheme. 

Algorithm: Adaptive quadrature construction 

Input: The domain of integration Q; numf integrands {/*}; prescribed error tolerance tol 
Output: Adaptive quadrature over Q that integrates {fi}™=™^ within the accuracy of tol 
for i — 1 : numf do 

11 = integrate fa over Q with 5-point quadrature 

1 2 = integrate fa over Q with 8-point quadrature 

if |/2 — h | > tol then set Pi — 1, otherwise set Pi = 

end do 

if {P l )ZT f = o 

return the 5-point quadrature as the rule over Q 

else 

partition Q into eight uniform cells 

for each cell Qj, call the quadrature construction routine with the integrands 

{/,•• S.t. = 

put the eight obtained quadratures together: the adaptive quadrature over Q 

end if 

This adaptive integration scheme is similar to that of van Dooren and de Ridder [28], 
except that in Reference [28], the domain is subdivided in one direction (dimension) at a time, 
and all integrands are carried until the last iteration. Limiting subdivision to one dimension 
at each step can be useful if the integrands have a strong dependence on one of the spatial 
dimensions. Berntsen et al. [29] devise an analogous scheme with nonuniform subdivision, 
and improve the error estimate: the error is approximated by combining integrals over a cell 
and its children. Adaptive integration by octree subdivision is also used by Pieper [30], and 
the integration error is estimated as the difference between the integral over a cell and its 
children. 

2.1 Example 

Consider the functions A(x) = 10exp(— 100r^) and ^(x) = 100 exp(— 200^) as the inte- 
grands for the quadrature construction, where T\ is the distance from the point x to the 
origin, r 2 is the distance to (.81, .62, .73), and the domain of integration is the unit cube. 
Figure la shows an adaptive quadrature over the unit cube with the above integrands and an 
error tolerance of 10 -6 with 8875 integration points. A sequence of octree refinements of the 
domain is shown in Figure lb-le. The MATLAB™ implementation of the adaptive quadrature 
is listed in Appendix A. 
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(a) (b) (c) (d) (e) 

Figure 1: Adaptive quadrature over the unit cube, (a) quadrature points; and (b)-(e) the 
sequence of octree refinements. 

3 Error Analysis 

In this section, the convergence properties of the adaptive integration algorithm are studied. 
Specifically, we are interested in cases where the integrand has a cusp, for example /(x) = 
1— r or /(x) = exp(— ar), which are C° functions with a derivative-discontinuity at the origin. 
However, this will not have an adverse effect on the efficiency of the adaptive integration 
scheme, and a high rate of convergence is realized. The following examples clarify the 
problem. 

Consider the integration of the function /(x) = 1— r, where r is the distance of the point x 
to the origin. The integration domain is [— 1, l] n for n = 1, 2, . . . , 6. In Figure 2, the function 
/ is illustrated in one and two dimensions. Gauss quadratures are used to integrate /, and the 
position of the cusp is arbitrary. By contrast, one could consider the cusp as a hindrance to 
convergence, and partition the integration domain by placing the cusp at the boundary of the 
subdivisions. Figure 3 shows the convergence of the integration as the number of integration 
points is increased: a rate of convergence of n + 1 is observed in M n , which can be explained 
by appealing to the Taylor expansion of the integrand. Since the integrand is continuous 
over a small cell Qq containing the cusp, the error of approximating it with a polynomial 
is 0{h)\ hence the error in the integration is O(h)Vn ~ 0{h nJrl ) ) where Vq ~ 0{h n ) is 
the volume of Qq. This indicates that our quadrature construction scheme can attain higher 
convergence rates in higher dimensions (for C° functions), in contrast to the usual deficiency 
of numerical integration methods that suffer from the curse of dimensionality. We perform 
the same experiment with /(x) = exp(— 20r), which has a cusp at the origin. The function 
is plotted in Figure 4a in two dimensions. The convergence curves for the integration using 
tensor-product quadratures are shown in Figure 4b with respect to the minimum distance of 
the integration points to the cusp. As expected, a rate of convergence of n + 1 is observed 
in W 1 . Numerical experiments also reveal that for smoother functions C m , m > 0, higher 
convergence rates are achieved, consistent with Darboux's Principle [31]. 
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(a) (b) 
Figure 2: /(x) = 1 — r. (a) in one dimension; and (b) in two dimensions. 
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Figure 3: Convergence curves of the integration of /(x) = 1 — r. (a) with respect to the 
number of integration points in each direction; and (b) with respect to the minimum distance 
of the Gauss points to the cusp (only even number of integration points are shown). 
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(a) (b) 

Figure 4: /(x) = exp(— 20r). (a) / in two dimensions; and (b) convergence curves of the 
integration of /(x) with respect to the minimum distance of the Gauss points to the cusp 
(only even number of integration points are shown). 

4 Numerical Examples 

4.1 Integration of the Regularized Heaviside Function 

Solution fields with sharp spatial gradients arise in modeling physical phenomena such as 
shear band evolution, damage and cohesive process zones, and convection-dominated prob- 
lems with shocks. Different classes of the regularized step function are adopted as the enrich- 
ment function for modeling sharp gradients, among which the following piecewise polynomial 
regularized Heaviside function is proposed by Patzak and Jirasek [1]: 

r if (j) < -e 

#M= ^/-,(l-S)V if |0|<e (1) 
[ 1 if (ft > e, 

where 0(x) is the signed distance from the interface, and e is a parameter that determines 
the gradient of the enrichment function (half- width of the zone) . The reference volume V £ is 
set to 256£:/315 to enforce C 4 continuity. Integrating (1) gives [10]: 

M, e) = — ^ (128s 9 + 315(f>e s - 4200V + 3780V - 1800V + 350 9 ) , (2) 
25o£ y 

in the region |0| < e. It has been suggested that a single enrichment function may not be 
adequate to represent the complete range of sharp gradients present in a problem, e.g., the 
shock front in a convection-dominated problem, or the damage process zone in quasibrittle 
materials. For such problems, multiple enrichment functions, each having a separate pa- 
rameter e should be used [1,10]. For the sake of illustration consider five functions with 
the parameter e/h = {2.5,0.85,0.265,0.085,0.0225} (h is the element size), and integrate 
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them on the unit square (h — 1). This set was produced by Abbas et al. [10] to minimize 
the point wise error in modeling shocks. Other types of sharp-gradient enrichment functions, 
such as tanh(g0) can also be used, where q controls the severity of the gradient. The interface 
is a straight line shown in Figure 5a, and the enrichment function for e = 0.085 is depicted 
in Figure 5f. Four different integration strategies are examined: (1) tensor-product over the 
square, ignoring the interface; (2) and (3) integration over subtriangles using a coarse and 
a fine partitioning; and (4) adaptive quadrature. In all cases, the accuracy is improved by 
increasing the number of integration points (Figure 5g). The adaptive integration outper- 
forms all other integration strategies. The experiment is repeated for a kinked interface and 
the results are shown in Figure 6. Due to the shape of the interface, the integrand in this 
case is relatively more complicated, and the advantage of using the adaptive quadrature is 
emphasized. For an arbitrary curved discontinuity, different methods have been proposed, 
for example, integration by partitioning the element into triangles with curved edges [32, 33], 
and octree-subdivision [34]. In Figure 7, the regularized Heaviside functions are integrated 
over a domain with a curved discontinuity. The interface is represented in closed form us- 
ing a quadratic polynomial. The performance of the adaptive quadrature is compared with 
the tensor-product quadrature in Figure 7e. Adaptive quadratures require fewer integration 
points for the same accuracy. 

4.2 Enriched Finite Element in Quantum- Mechanical Calculations 

In Reference [14], an enriched finite element method was applied to the Coulomb problem 
in crystalline diamond, and tensor-product quadratures were used for the numerical integra- 
tion. The number of integration points was increased until convergence in the solution was 
achieved. A higher-order quadrature was used over the elements containing the nuclei, since 
the enriched basis functions and the source term were strongly localized about the nuclei. 
A large number of integration points were required to obtain the desired accuracy and the 
optimal rate of convergence. Herein, we use the adaptive integration scheme introduced 
in Section 2 to setup the system matrices, and show that the optimal rate of convergence 
is realized at a relatively low computational cost. First, we present a brief description of 
the problem, and then the adaptive numerical integration algorithm is used to set up the 
finite element system matrices. For details on the formulation and solution technique, see 
References [13, 14]. 

Consider the unit cell defined by the lattice vectors 



where a = 6.75 bohr, and carbon atoms are at T\ = (0,0,0) and r 2 = (1/4,1/4,1/4) in 
lattice coordinates. The total charge p in the unit cell is written as 



ai = |(0,l,l) 
a 2 = |(1,0, 1) 

83 = 1(1,1,0) 



(3) 



p(x) = p+(x) + p-(x) = p + (x) - p(x) + p-(x) + p(x) = p+(x) + p-(x), 



(4) 
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(a) (b) (c) (d) (e) 
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(f) (g) 

Figure 5: Numerical integration of the regularized Heaviside function with a straight inter- 
face, (a) domain of integration and interface; (b) tensor-product; (c) triangular quadratures 
(4 triangles); (d) triangular quadratures (8 triangles); (e) adaptive quadrature; (f) the regu- 
larized Heaviside function for e = 0.085; and (g) maximum relative error in the integration 
of the five functions versus the number of integration points for different strategies. 
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(a) (b) (c) (d) (e) 




Total number of integration points 



(f) (g) 

Figure 6: Numerical integration of the regularized Heaviside function with a kinked inter- 
face, (a) domain of integration and interface; (b) tensor-product; (c) triangular quadratures 
(6 triangles); (d) triangular quadratures (32 triangles); (e) adaptive quadrature; (f) the reg- 
ularized Heaviside function for e = 0.085; and (g) maximum relative error in the integration 
of the five functions versus the number of integration points for different strategies. 
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Total number of integration points 



(d) (e) 

Figure 7: Numerical integration of the regularized Heaviside function with a curved inter- 
face, (a) domain of integration and interface; (b) tensor-product; (e) adaptive quadrature; 
(f) the regularized Heaviside function for e = 0.085; and (g) maximum relative error in 
the integration of the five functions versus the number of integration points for different 
strategies. 
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Figure 8: (a) Electronic charge densities; and (b) enrichment functions for the Coulomb 
problem in crystalline diamond [14]. 



where p + (x) = ^ Pi(x) = %5(x — Tj) is the total nuclear charge density in the unit cell, 
p~(x) is the electronic charge density, and p + (x) = p + (x)— p(x) and p~(x) = p~(x)+p(x) are 
the neutralized nuclear and electronic charge densities, respectively. The neutralizing charge 
p(x) is introduced in (4) to circumvent the divergence of the potential ^ + (x) (V + ~ 1/r) at 
nuclear locations, so that V + is extracted analytically, and V~ is solved in real space. The 
potential V~ associated with the neutralized electronic charge density is obtained from the 
solution of Poisson's equation: 

V 2 y-(x) = -47rp-(x), (5) 

subject to periodic boundary conditions, with continuous neutralized electronic charge den- 
sity p~(x) as the source term. The EFE solution is written as 

y-(x) = ^<&(x)ai + $^a(x)6 a = ^$ fc (x)c fe , (6) 

i a k 

where a is summed over the atoms and = {4>i} U {tp a } is the combination of the 

classical and enriched basis functions that form the EFE basis. The enrichment functions 
^a(x) are taken as sum of the potentials vj — isolated atomic solutions corresponding to the 
neutralized electronic charge densities pj = pj + pj in the vicinity of each atom /. The 
enrichment function is written as 

^(x) = y-(x) = J>-(x-R), (7) 

R 

where R denotes lattice translation vectors. The electronic charge densities and the enrich- 
ment functions are plotted in Figure 8. 

On incorporating the trial and test functions of the form (5) into the weak form of the 
Poisson's equation, the discrete linear system of equations emerges. It is seen that the terms 
ijj,ip 2 ,dilj/dr,(dip/dr) 2 , and ipdi/j/dr appear in the element stiffness matrix, and p in the 
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Figure 9: Normalized integrands in the element stiffness matrix and force vector of the 
Poisson problem. 



element force vector, where for simplicity ip and p are used for the enrichment function, and 
the electronic charge density, respectively, and dip/dr refers to the derivative with respect 
to the radial coordinate. These functions are normalized, with an absolute maximum value 
of unity, and plotted in Figure 9. The integrands have sharp gradients close to the atomic 
positions and produce cusps at the atomic sites. While it is possible to integrate these 
terms separately (resulting in multiple ad hoc quadratures), it is desirable to have a single 
quadrature that is capable of efficiently evaluating the integrals altogether. The proposed 
numerical integration algorithm constructs a quadrature rule over each finite element that 
satisfies a given error tolerance for all the above integrands. Once the finite element mesh is 
generated, adaptive quadratures are constructed and saved for each element, which are used 
throughout the analysis and the postprocessing. 

A sequence of refined meshes are used with 4, 8, 12, 16, 20, 24, and 32 cubic serendipity 
finite elements in each direction. Numerical integration is performed with tensor-product 
and adaptive quadratures. In each case, the accuracy of the quadratures is increased until 
convergence in the solution is observed. The number of integration points for each mesh 
and integration strategy is reported in Table 1. Figure 10 shows the convergence curves for 
FE and EFE runs, with a rate of convergence of 4.39 and 5.96, respectively, for the last 
three data points. The results indicate that the integration scheme is sufficiently accurate 
to realize the optimal rate of convergence. The error tolerance (input of the quadrature 
construction algorithm) is the maximum allowable error that produces a stable result (i.e., 
stable with respect to further decrease in tolerance). For the pure FE runs, a 5-point Gauss 
quadrature rule is used in each direction. Adaptive integration proves to be superior with 
respect to the tensor-product quadrature. The improvement is emphasized for finer meshes 
where higher accuracies are required: the integration demand of the EFE solution is only 
marginally higher than the pure FE solution (of much lower accuracy) on the same mesh. 
In Figure 11, a comparison is made between the behavior of tensor-product and adaptive 
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Table 1: Comparison of the number of integration points for tensor-product and adaptive 
quadratures. 



Mesh 


Error tolerance 


Number of integration points 








Pure FE 


Tensor-product 


Adaptive 


4 


2.2 x 10" 3 


8000 


169000 


78000 


8 


1.7 x 10" 4 


64000 


624000 


162000 


12 


2.2 x 10" 5 


216000 


1840000 


349000 


16 


4.3 x 10" 6 


512000 


14020000 


981000 


20 


1.1 x 10" 6 


1000000 


27196000 


1549500 


24 


3.0 x 10" 7 


1728000 


46852000 


3073750 


32 


4.7 x 10" 8 


4096000 




6112000 



quadratures for the EFE solution of the Coulomb problem in crystalline diamond. The 
quadrature error in the Coulomb energy is plotted with respect to the number of integration 
points for finite element meshes with 4 and 16 elements in each direction. The accuracy 
of the quadrature is increased until the desired convergence with respect to the quadrature 
is attained. The adaptive quadrature requires fewer integration points to achieve the same 
accuracy. Standard Gauss quadrature shows a smoother convergence curve, which can be 
attributed to the uniform overall increase in the accuracy of the quadrature throughout 
the domain. The efficiency of the adaptive quadrature is more noticeable in case of the 
16 x 16 x 16 mesh. 

Note that in all the meshes used earlier, the atoms are located at the vertices of the 
finite elements — the atom inside the unit cell is at a quarter of the diagonal from the corner, 
and the number of elements is a multiple of four — and there is no atom inside any finite 
element. For more general lattice systems, non-uniformly refined meshes may be required in 
order to ensure that all atoms are located at the vertices of elements. This is not desirable 
due to the increase in the total number of elements and the associated computational costs. 
Furthermore, with the application of the EFE, one would like to resolve the local features 
by adding enrichment functions to the approximation, and not by refining the finite element 
mesh. In the following, we use meshes with 1, 2, 3, 5, 6, and 7 elements in each direction. In 
all these cases, the atom in the unit cell lies inside an element. The number of integration 
points of the tensor-product and adaptive quadratures are compared in Table 2. The increase 
in the number of integration points for the tensor-product quadrature is significant, whereas 
the adaptive integration provides accurate results with a moderate number of integration 
points. 

5 Concluding Remarks 

An adaptive integration scheme was presented that can be used for functions with sharp 
gradients and cusps. The algorithm uses tensor-product quadratures with 5 and 8 integra- 
tion points in each direction to estimate the local error, and divides the domain uniformly, 
independent of the location of the cusp or sharp gradient, until the target tolerance is met. 
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Number of elements in each direction 



Figure 10: Error in Coulomb energy per atom. 32 x 32 x 32 mesh is used as the reference 
solution. 




,8 



Number of integration points Number of integration points 



(a) (b) 
Figure 11: Comparison of tensor-product and adaptive quadratures for the crystalline dia- 
mond problem. The thick line shows the EFE solution error, (a) 4 x 4 x 4 mesh; and (b) 
16 x 16 x 16 mesh. 



15 



Table 2: Comparison of the number of integration points for tensor-product and adaptive 
quadratures: one of the atoms is inside an element. 



Mesh 


Error tolerance 


Number of integration points 








Pure FE 


Tensor-product 


Adaptive 


1 


9.8 x 10" 3 


125 


1000000 


48250 


2 


6.8 x 10" 3 


1000 


4096000 


48250 


3 


3.3 x 10" 3 


3375 


1880000 


55875 


5 


1.1 x 10" 3 


15625 


8000000 


165250 


6 


4.9 x 10" 4 


27000 


27000000 


188875 


7 


3.0 x 10" 4 


42875 


21952000 


183750 



The error analysis in the integration of a function with a cusp (derivative-discontinuity at 
a point) showed that the rate of convergence is improved in higher dimensions. The adap- 
tive integration algorithm was successfully used for the integration of a set of regularized 
Heaviside functions, and proved to be more efficient than integration by partitioning as well 
as tensor-product quadratures. The method was also applied to the enriched finite element 
solution of the all-electron Coulomb potential and energy of crystalline diamond (Poisson's 
equation). The enrichment functions and source term were strongly localized about the 
atomic positions. The adaptive integration scheme proved to be very efficient, and recovered 
the optimal rate of convergence with only a moderate increase in the number of integration 
points with respect to the classical finite element method (of much lower accuracy) on the 
same mesh; while reducing the integration points required in the EFE solution by an order 
of magnitude or more. The adaptive integration scheme is simple, robust, and directly ap- 
plicable to any generalized finite element method employing enrichments with sharp local 
variations or cusps in n-dimensional parallelepiped elements. 
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A MATLAB Code for the Adaptive Integration Scheme 

The following routine is the implementation of our adaptive integration scheme and produces 
a quadrature over an n-dimensional hyperparallelepiped for a given set of functions and 
prescribed error tolerance. First, a description of the input and output parameters is given. 

• fn: 1 x numf , set of integrands, defined as an array of structures with the mem- 
ber h, which is a function handle. For example, in case of two integrands, we have: 
fn(l).h = ©integrandl ; and fn(2).h = @integrand2; , where integrandl .m and 
integrand2.m are MATLAB functions defined as W 1 —> R. 

• d: (n + 1) x n, domain of integration, a hyperparallelepiped in W 1 , defined using one 
point as its base and only n subsequent vertices of d. For example, while d(2, : ) is a 
vertex of d, the vector d(2, : )-d(l , : ) is a lattice vector of d. 

• nsp: 1x2, number of integration points in each direction to evaluate the integral over 
the partitions, nsp(l) is used to evaluate the local integrals, and nsp (2) is used to 
evaluate the local integration error. 

• tol: lxl, absolute value of the integration error. 

• X: n x numx, quadrature points, each point is a column vector. 

• W: numx x 1, quadrature weights. 

For example, the functions used in Section 2.1 are defined in the following m-files: 

function val = integrandl (X) 
R = sqrt(sum(X. ~2, 1)) ; 
val = 10*exp(-100*R. ~2) ; 
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% 

function val = integrand2(X) 

R = sqrt(sum((X - repmat ( [. 81 , .62, .73]', 1, size(X, 2))).~2, D); 
val = 100*exp(-200*R. ~2) ; 

The domain of integration (unit cube) is defined using its base and three of its vertices: 

d = [0, 0, 0; . . . 
1, 0, 0; ... 
0, 1, 0; ... 
0, 0, 1]; 

With the above setting, the adaptive quadrature function can be called (see Figure 1): 
[X, W] = ndimensional_adaptive_integration(f n, d, [5, 8], le-6) ; 
The MATLAB code for the quadrature construction follows. 

function [X, W] = ndimensional_adaptive_integration(f n, d, nsp, tol) 

% External Dependencies (m-files) 

% gauss_points(nsp) : ID Gauss points 

% gauss_weights(nsp) : ID Gauss weights 

dim = size(d, 2); % dimension 

[XYZgl, Wgl] = Get_initial_quad(nsp(l) , dim); 

[XYZg2, Wg2] = Get_initial_quad(nsp(2) , dim); 

[X, W] = Adaptive_integration(f n, d, dim, . . . 

nsp, XYZgl, Wgl, XYZg2, Wg2, tol); 
% 

function [XYZ, W] = Adapt ive_integrat ion (f n, d, dim, . . . 

nsp, XYZgl, Wgl, XYZg2, Wg2, tol) 

numfun = length (fn) ; 

nspd = nsp. "dim; 

PARTITION = zeros(numfun, 1); 

for ind = 1: numfun 

[integl, XYZtl, Wtl] = Integrate_over_one_cell(fn(ind) , ... 

d, dim, nspd(l), XYZgl, Wgl); 
integ2 = Integrate_over_one_cell(fn(ind) , d, dim, ... 

nspd(2), XYZg2, Wg2) ; 

err = abs(integ2 - integl); 

if (err >= tol), PARTITION (ind) = 1; end 

end 

if any (PARTITION) 

two_power_dim = 2~dim; 

XYZtemp = cell(l, two_power_dim) ; Wtemp = XYZtemp; 
vectors = (d(2:end, :) - repmat (d(l, :), dim, 1)) / 2; 
I = zeros(dim, 1) ; 
for ind = 1 : two_power_dim 
base = d(l , : ) ; 
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for j = l:dim 

base = base + (I(j))*vectors(j , :); 

end 

division = [base; repmat (base , dim, 1) + vectors] ; 
[XYZtemp{ind}, Wtemp{ind}] = Adapt ive_integrat ion ( . . . 

fn(PARTITION == 1), division, dim, nsp, ... 

XYZgl, Wgl, XYZg2, Wg2, tol) ; 

I = incrementl (I , dim); 

end 

XYZ = [] ; W = [] ; 

for ind = 1 : two_power_dim 

XYZ = [XYZ, XYZtemp{ind}] ; 

W = [W; Wtemp-Cind}] ; 

end 

else 

XYZ = XYZtl; 
W = Wtl; 

end 

% 

function [integ, XYZt, Wt] = Integrate_over_one_cell(fn, d, ... 

dim, nspd, XYZg, Wg) 
A = (d(2:end, :) - repmat (d(l, :), dim, 1)) ; ; 
shift = d(l, :)>; 

XYZt = A*XYZg + repmat (shift, 1, nspd); 

Wt = Wg*det(A) ; 

integ = dot(Wt, fn.h(XYZt)); 

% 

function [XYZg, Wg] = Get_initial_quad(nsp , dim) 

% Get a nsp~dim point quad over the unit 'dim' -dimensional hypercube 

xg = gauss_points(nsp) ' ; xg = (1 + xg)/2; 

wg = gauss_weights(nsp) ' ; wg = wg/2; 

XYZg = xg; 

Wg = wg; 

for i = 2: dim 

XYZg = repmat (XYZg, nsp, 1); 

XYZg(:, i) = reshape (repmat (xg, 1, nsp~(i-l)) ; , 1, nsp~i) ; 
Wg = repmat (Wg, nsp, 1); 

Wg(:, i) = reshape (repmat (wg, 1, nsp~(i-l)) ; , 1, nsp~i) ; 

end 

XYZg = XYZg' ; 

Wg = prod(Wg, 2) ; 

% 

function I = incrementl (I , dim) 

% Keep track of the vectors that construct the current cell 
1(1) = 1(1) + 1; 
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for i = 1 : (dim-1) 

if (I(i) == 1), break; end 
Ki) = 0; 

I(i+1) = I(i+1) + 1; 

end 



